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Abstract. We perform a careful analysis of the main Monte Carlo algorithm used in parton shower sim- 
ulations, the Sudakov veto algorithm. We prove a general version of the algorithm, directly including the 
dependence on the infrared cutoff. Taking this as a starting point, we then consider non-positive definite 
splitting kernels, as encountered when dealing with sub-leading colour correlations or splitting kernels 
beyond leading order. New algorithms suited for these situations are developed. 

PACS. 02.70.Tt Monte Carlo methods - 12.38.Cy Summation of QCD perturbation theory 
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1 Introduction 

Parton shower Monte Carlo simulations as implemented 
in for example [IH3], are indispensable tools for analyzing 
and predicting realistic final states encountered in collider 
experiments. Matrix element corrections, as discussed in 
[3HH]; the technically similar matching to NLO calcula- 
tions employing the POWHEG method [9], or schemes to 
combine parton showers and multijet tree-level matrix el- 
ements |10H14j . all rely directly or indirectly on the same 
method for generating subsequent parton shower emis- 
sions in Monte Carlo simulations. 

With the notable exception of the FORTRAN version 
of HERWIG, nowadays most parton shower implementa- 
tions use the Sudakov veto algorithm to facilitate this task, 
as the splitting kernels normally are too complicated to al- 
low efficient integration. 

A justification of the Sudakov veto algorithm is given 
in [1] , stating that for upper bounds R on splitting kernels 
P, R(q) > P(q) for all q, algorithm (JTJ) will draw random 
variables with density 

dT P (q\Q) = 0(Q - q)P(q)A P (q\Q)dq , (1) 

where the Sudakov form factor is given by 

A P {q\Q) = cxp (- £ P®**) ■ (2) 

We note here, however, that the algorithm has to be 
more carefully defined. Most obviously if, in algorithm (JXJ) , 
P(q) = but R(q) ^ for all q < q c and some q c , the 
algorithm will potentially enter an infinite loop. We shall 
therefore assume that R(q) is suitably restricted to avoid 
this situation, making the algorithm well-defined in the 
sense that it will never hit a state in which it will not 
terminate with probability one. 



Algorithm 1 The Sudakov veto algorithm as quoted in 
the literature. 



Q' 

loop 

Draw q with density 



9{Q' - q)R(q)A R (q\Q')dq 



return q with probability P(q)/R(q) 
Q'^q 
end loop 



Literally implementing the algorithm as presented above 
will not generate the desired density owing to the fact that 
dTp is not a probability density, 



Q dT P (t\Q) 



dt 



dt = 1- A P (q\Q) + 1 



(3) 



At best the algorithm will approximate the target den- 
sity if, for the lowest possible q, Ap{q\Q) «C 1. In practice, 
however, a vanishing Ap will never be encountered in par- 
ton shower simulations, due to the fact that an infrared 
cutoff [i > is always present. Thus the typically diver- 
gent part of the splitting kernel at q = is never reached, 
and the no-emission probability remains, Ap(fi\Q) > 0. 

Similarly, the competing processes algorithm 



Draw {qi, q n } from dT Pi (qi\Q), i 
return max({gi, q„}) 



targeting at drawing random variables with density dTp, 
P = J2i-Pi> wm n °t produce the desired result for the 
same reason. 
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2 The Complete Algorithm 

The failure of the simple algorithm presented in the pre- 
vious section has been argued to originate from the fact 
that the density considered is not a probability density. 

However, the density considered in the previous section 
is also not what is typically aimed at in a parton shower 
implementation. (See e.g. |15) for a concise treatment.) 
This can be seen by the fact that a lower cutoff scale /x 
has not been specified, nor is a virtual no-emission contri- 
bution present. Owing to the fact that a parton shower is 
to preserve the total inclusive cross section, the combined 
density, including both emission and no-emission, has to 
be a probability density. As the probability to not emit 
between two scales is determined by the Sudakov form 
factor, the probability density which we are interested in 
is 



A P (p\Q)6{q - n) + 

- q)9(q - l x)P{q)A P {q\Q) , (4) 



dS P (»,g\Q) 
dq 

which relates to the previously introduced density as 
dS P (n,q\Q) _ 



dq 



A P (ji\Q)6(q -n) + 9(q - fj,)dT P (q\Q) . (5) 



Using sampling by inversion, Q 



q dS P (n,t\Q) 
dt 



dt = A P (q\Q)9(q - fj,) = rnd 



(6) 



we find the equation to be solved for the next scale q. 
This is similar to what one would expect by viewing the 
Sudakov form factor Ap(q\Q) as a no-emission probability 
between two scales Q and q, but now explicitly taking into 
account the dependence on the infrared cutoff /i. 

As the splitting kernel, P, is not normally easily inte- 
grated, what is used in actual implementations is instead 
typically a version of the Sudakov veto algorithm where 
A R (q\Q) = rnd is solved for some easily integrable func- 
tion R(q) > P{q), and the radiation is kept only with a 
probability of P(q) / R(q). The issue of how to deal with 
the fact that the Sudakov factor A R (fi\Q) =t 0, however, 
remains. 

In the typically encountered case that P(q) is divergent 
at an absolute lower bound (which we take to be q = 0), 
the problem with the non- vanishing Sudakov factor at the 
lowest physically considered bound (q = /i) can be circum- 
vented by integrating down to q = 0. Events which only 
have emissions below the lowest physical bound (p) are 
then regarded as no-emission events [T^lfTT] . This is guar- 
anteed to work as for such splitting kernels Ap(0\Q) = 0. 

However, for splittings of massive particles it is the 
case, that - even if the splitting kernel is integrated down 



1 In this paper rnd denotes a source of uniformly distributed 
random numbers on [0, 1). 



to - the corresponding Sudakov factor is not vanishing, 
Ap(0\Q) 7^ 0. This situation can be dealt with by using 
an overestimation function R(q) which does correspond 
to A R (0\Q) ~ 0. The approximation of a non-divergent 
splitting kernel with a divergent one is, however, likely to 
lead to a severe overestimate, i.e. R(q) ^> P(q), which sig- 
nificantly influences the efficiency of the algorithm. Alter- 
natively, we here suggest that algorithm @ can be used, 
both for divergent and non-divergent splitting kernels. 

Algorithm 2 The alternative Sudakov veto algorithm. 
Q' 

loop 

solve rnd= Au(q\Q')8(q — /i) for q 
if q = fi then 

return fi 
else 

return q with probability P(q)/R(q) 
end if 
Q' ^q 
end loop 



We claim that this algorithm will correctly produce 
dSp([i, q\Q) for all chosen boundaries [i < q < Q. To 
prove it, we first prove theorem ([T]). 

Theorem 1 The q-density produced by the Sudakov veto 
algorithm after n rejection steps and a final acceptance 
step is given by 



dSi:l(^q\Q) 



A R ^\Q)5(q-^A { ;l R (^\Q)+ 



dq 



- q)9(q - t,)P(q)A R (q\Q)A^_ R (q\Q) (7) 



where 



a { pUi\Q) = ^ 



(P(k) - R(k))dk 



(8) 



From this the correctness of the algorithm follows upon 
summing over any number of rejection steps n = to oo, 
and the usage of A R (q\Q)A P - R (q\Q) = A P {q\Q). Note 
that theorem ([T]) does include the density of non-radiating 
events, and that each time the loop in algorithm is 
entered, an event q is drawn from dS R by virtue of eq. (J6]). 

We will show theorem (JT|) using induction and there- 
fore start by noting that the probability to accept an 
event, starting at an intermediate scale fc, is given by 

dS^ t (^q\k) = A R ( f ,\k)d(q-^)+ 

9(k-q)9(q- f x)P(q)A M (q\k) , (9) 

where the first term reflects the fact that proposal events 
at the infrared cutoff are always accepted, while the sec- 
ond term accounts for proposal events above the cutoff 
being accepted with probability P(q)/R(q). For n = the 
intermediate scale k equals the starting scale Q, i.e. 



(10) 
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proving eq. Q for n = 0. If the algorithm had performed 
one rejection step, events could only have been proposed 
above the infrared cutoff (otherwise the algorithm would 
have terminated), and we have 



dS£to(M)?|Q) 



dS acccp V q\k) (R(k) - P{k)) A R {k\Q)dk 
d41( M) #) ( R ( k ) - p ( fc )) A R (k\Q)dk . 



ill) 



Here, the factor of R(k) — P(k) originates from the veto 
probability, 1 — P(q)/R(q), times the kernel R(q) which 
had been used for the proposed event. 

To arrive at the desired density in eq. ((7} we use the 
'chain' property of the Sudakov form factors, 

A R (q\k)A R {k\Q) = A R {q\Q) , (12) 

and the relation 

rQ 

A { pl R {q\Q) = / (R(k) - P(k)) dk . (13) 

Jq 

This proves eq. ([7J for n = 1. In general, 

dS%+ 1) (n,q\Q) = 

&S%L(n, q\k) (R(k) - P(k)) A R (k\Q)dk (14) 

reflecting an initially proposed event k below Q, which 
initiated a sequence of n veto steps and a final acceptance 
step. Thus, if the theorem was correct for some n > 0, we 
readily obtain the claimed result for n + 1 upon using 




f(k')dk' f(k)dk 



ri+l 



(rc+1)! 



/(fc)dfc 



(15) 



The competing processes algorithm in turn reads 



Draw {qi, q n } from dSp^q^Q), i = 1, n 
return max({gi, q n }) 

which is easily proven as dSp i (qi\Q) now is a true proba- 
bility density. 



3 Towards Splitting Kernels of Indefinite Sign 

For the remainder of this note we shall be concerned with 
seeking solutions to the case of non-positive definite split 



or D(x) < 0, however, the generalization of the Sudakov 
veto algorithm is not obvious. 

In this section we will outline an algorithm, algorithm ([3]) , 
which is able to deal with the general case of non-positive 
definite splitting kernels, but is limited to considering dis- 
tributions at fixed starting scale Q. That this is a limita- 
tion can be seen from the fact the the generated density 
will multiply a Q-dependent normalization factor smaller 
than one. As long as only one starting scale is considered, 
this scale dependence can trivially be normalized away. 
However, in the case of varying scales, one would have to 
introduce scale dependent event weights larger than one. 
For the case of an unlimited number of emissions driven by 
subsequently sampling the density at varying scales, there 
is clearly no upper bound for the combined size of these 
weights. The algorithm presented here could, however, be 
of practical interest for cases where splitting kernels P 
of indefinite sign are present only for a limited number 
of emissions. Such scenarios would indeed give rise to an 
upper bound on the expected event weight; particularly 
one could consider matrix element corrections incorporat- 
ing higher order corrections with the need for appropriate 
subtractions to regularize infrared divergences. 

To be precise, we decompose the non-positive definite 
kernel P(q) as P + (q) — P~(q), where 



P 



±()= (±P(q) : P(q)^0 
I : otherwise 



(16) 



and utilize algorithm ([3]). 



Algorithm 3 The algorithm for splitting kernels of indef- 
inite sign. See text for the definition of P±. 
loop 

Draw q+ from dS P + (/i, q\Q) 
Draw q- from dS P - (fi, q\Q) 
q <- max(q+,q-) 
if q = fi then 

return /i with weight +1 
end if 

Draw t from dS 2 p- {n, t\q) 
if t — n then 

if max(g+ , <?- ) = q+ then 
return q with weight +1 

else 

return q with weight —1 
end if 
end if 
end loop 



Note that all random variables needed from Sudakov- 
type distributions are readily generated using the veto al- 



ting kernels. For potentially negative- valued 'densities' J D(x),g° rithm as outlined above. We claim that the algorithm 
a Monte Carlo implementation is still sensible by sam- w ^ g enera te 
pling events x according to \D{x)\ and afterwards assign- 
ing weights +1 or —1, depending on whether D(x) > dSp(fx, q\Q) x Ap-([i\Q) . 
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To prove this, note that if q+ = q- — /i we obtain a 
contribution 

S{q-n)A P +Ot\Q)Ap-(ji\Q) = 

5{q- ii)A P {n\Q)A%-(n\Q) . 

The probability for t — fi, i.e. not to re-enter the loop is 
clearly given by Ap_(fi\q). Then, if q + > g_ (and hence 
q+ > fi), we hnd a contribution 

6{q - /i)dT P+ (q\Q)A P - (q\Q)A 2 P - (jj,\q) = 

9(q - M )dT P+ (q\Q)A_ P - (q\Q)A 2 p - (fi\Q) . 

Finally, if g_ > q+ (and hence g_ > /x), while including 
the negative weight for these events, the last contribution 
is 

- % - fi)dT P - (q\Q)A P+ (q\Q)A 2 P - (ji\q) = 

9(q - v)dT_ P - (q\Q)A P+ (q\Q)A 2 p - (fi\Q) , 

completing the proof. For the case of several available pro- 
cesses we can always decompose 

p(q) = J2 p m = E^ + (?) -J2 p r(i) > ( 17 ) 

i i i 

such that the proposal events q±, as well as the 'control 
variate' t may be generated using the competing processes 
algorithm for the individual positive and negative contri- 
butions, 

P ± {q)=J2 P t(l)- (18) 



4 Interleaving Vetoing and Competition 

The algorithm outlined in the previous section may be 
used to deal with the case of non-positive definite split- 
ting kernels in full generality provided we are interested 
in distributions for a single starting scale Q, or are pre- 
pared to accept potentially large weights. For practical 
purposes, we are, however, interested in cascades at sub- 
sequent scales qi > q2 > ■■■ > q n , where qk-i serves as the 
starting scale of the distribution for qf. . The Q-dependent 
normalization Ap_ present in the distribution generated 
will thus make it non-ideal in the context of cascades. 

Here we consider the typically encountered physical 
setup for which we may assume that 

P(q)=J2 P M>0, (19) 

i 

still allowing for a probabilistic interpretation, though a 
subset of the splitting kernels are of indefinite sign. P(q) 
can be decomposed as in eq. (jTTJ) , and we can directly 
identify an overestimate to the desired splitting kernel, 

P + (q)>P(q)=P + (q)-P-(q). (20) 



Algorithm 4 The interleaved veto/competition algo- 
rithm. 

Q' 

loop 

Draw {qi, ...,<?„} from dS p + (qi\Q), i = 1, —,n 

q <- max({qi,...,q n }) 
if q — fi then 

return /x 
else 

return q with probability {P + {q) — P~ (q)) / P + (q) 
end if 
Q' ^q 
end loop 



This suggests a two-step procedure of interleaving com- 
peting processes and vetoing, formalized in algorithm (j4|), 
which we choose to call the 'interleaved veto/competition 
algorithm'. Here, the qi may be generated directly, if the 
P^~ allow to. Alternatively the veto algorithm may be used 
with overestimates Rf(q) > P^~(q). The correctness of the 
complete algorithm is seen by the fact that the first two in- 
structions in the loop will guarantee that q is distributed 
according to dSp+ by the competing process algorithm. 
In the following steps, the obtained density P + {q) is cor- 
rected to P(q) = P + (q) — P~(q) by virtue of the standard 
veto algorithm. Note that this algorithm will neither re- 
quire negative weights, or introduce a Q-dependent nor- 
malization. 



5 Conclusions and Outlook 

We have given a careful analysis of the main Monte Carlo 
algorithm entering current parton shower simulations, the 
Sudakov veto algorithm. Especially, we have discussed in 
detail the importance of the no emission probability aris- 
ing as a consequence of an infrared cutoff, and suggested 
an alternative formulation, algorithm which directly 
includes the dependence on the infrared cutoff. This al- 
gorithm is argued to be more efficient in the case of a 
non-divergent splitting kernel. 

We also consider possible extensions to the case of 
splitting kernels of indefinite sign. Such splitting kernels 
are encountered when trying to extend parton showers be- 
yond the large N c limit or beyond leading order. 

First, in algorithm ([3]) we develop a general algorithm 
for a splitting kernel of indefinite sign. Modulo a nor- 
malization dependence on the starting scale of the algo- 
rithm, this case may indeed be dealt with in full gener- 
ality. The Q-dependent normalization, however, prevents 
efficient usage in the context of cascades using an ordered 
chain of scales. 

For the typically encountered case, in which splitting 
kernels of indefinite sign are present, but the sum over all 
possible splitting kernels stays positive, we give, in algo- 
rithm @ , an algorithm interleaving the competing process 
algorithm with subsequent veto steps. 
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